{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "验证 f_x 结果： 0\n",
      "验证 f_y 结果： 0\n"
     ]
    }
   ],
   "source": [
    "import sympy as sp\n",
    "\n",
    "# 定义符号变量\n",
    "x, y, t = sp.symbols('x y t')\n",
    "nu = 1  # 已知 ν=1\n",
    "\n",
    "# 速度分量、压力以及已知的 f\n",
    "ux = (x**2 * y**2 + sp.exp(-y)) * sp.cos(2 * sp.pi * t)\n",
    "uy = (-2/3 * x * y**3 + 2 - sp.pi * sp.sin(sp.pi * x)) * sp.cos(2 * sp.pi * t)\n",
    "p = - (2 - sp.pi * sp.sin(sp.pi * x)) * \\\n",
    "    sp.cos(2 * sp.pi * y) * sp.cos(2 * sp.pi * t)\n",
    "# 修正 y^2 为 y**2\n",
    "fx_given = -2 * sp.pi * (x**2 * y**2 + sp.exp(-y)) * sp.sin(2 * sp.pi * t) + (-2 * nu * x**2 - 2 * nu *\n",
    "                                                                              y**2 - nu * sp.exp(-y) + sp.pi**2 * sp.cos(sp.pi * x) * sp.cos(2 * sp.pi * y)) * sp.cos(2 * sp.pi * t)\n",
    "fy_given = -2 * sp.pi * (-2/3 * x * y**3 + 2 - sp.pi * sp.sin(sp.pi * x)) * sp.sin(2 * sp.pi * t) + (4 * nu * \\\n",
    "    x * y - nu * sp.pi**3 * sp.sin(sp.pi * x) + 2 * sp.pi * \\\n",
    "    (2 - sp.pi * sp.sin(sp.pi * x)) * \\\n",
    "    sp.sin(2 * sp.pi * y)) * sp.cos(2 * sp.pi * t)\n",
    "    # 这里缺cos(2*pi*t)\n",
    "\n",
    "# 计算速度对时间的导数 u_t\n",
    "ut_x = sp.diff(ux, t)\n",
    "ut_y = sp.diff(uy, t)\n",
    "\n",
    "# 计算速度梯度 ∇u\n",
    "d_ux_x = sp.diff(ux, x)\n",
    "d_ux_y = sp.diff(ux, y)\n",
    "d_uy_x = sp.diff(uy, x)\n",
    "d_uy_y = sp.diff(uy, y)\n",
    "nabla_u = sp.Matrix([[d_ux_x, d_ux_y], [d_uy_x, d_uy_y]])\n",
    "\n",
    "# 计算变形率张量 D(u)\n",
    "D_u = 0.5 * (nabla_u + nabla_u.T)\n",
    "\n",
    "# 构造应力张量 T = 2νD(u) - pI\n",
    "T = 2 * nu * D_u - p * sp.eye(2)\n",
    "\n",
    "# 计算 ∇·T\n",
    "div_T_x = sp.diff(T[0, 0], x) + sp.diff(T[1, 0], y)\n",
    "div_T_y = sp.diff(T[0, 1], x) + sp.diff(T[1, 1], y)\n",
    "\n",
    "# 计算理论 f（ut - ∇·T）\n",
    "fx_calc = ut_x - div_T_x\n",
    "fy_calc = ut_y - div_T_y\n",
    "\n",
    "# 验证给定的 f 是否正确（计算给定 f 与理论 f 的差）\n",
    "verify_fx = fx_given - fx_calc\n",
    "verify_fy = fy_given - fy_calc\n",
    "\n",
    "print(\"验证 f_x 结果：\", sp.simplify(verify_fx))\n",
    "print(\"验证 f_y 结果：\", sp.simplify(verify_fy))\n",
    "\n",
    "# sp.latex(fy_calc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - 2.0 x y^{2} \\left(- 0.666666666666667 x y^{3} - \\pi \\sin{\\left(\\pi x \\right)} + 2\\right) \\cos^{2}{\\left(2 \\pi t \\right)} + \\left(- 0.666666666666667 y^{3} - \\pi^{2} \\cos{\\left(\\pi x \\right)}\\right) \\left(x^{2} y^{2} + e^{- y}\\right) \\cos^{2}{\\left(2 \\pi t \\right)}$"
      ],
      "text/plain": [
       "-2.0*x*y**2*(-0.666666666666667*x*y**3 - pi*sin(pi*x) + 2)*cos(2*pi*t)**2 + (-0.666666666666667*y**3 - pi**2*cos(pi*x))*(x**2*y**2 + exp(-y))*cos(2*pi*t)**2"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 带时间ns\n",
    "import sympy as sp\n",
    "\n",
    "# 定义符号变量\n",
    "x, y, t = sp.symbols('x y t')\n",
    "\n",
    "# 定义 ux 和 uy\n",
    "ux = (x**2 * y**2 + sp.exp(-y)) * sp.cos(2 * sp.pi * t)\n",
    "uy = (-2/3 * x * y**3 + 2 - sp.pi * sp.sin(sp.pi * x)) * sp.cos(2 * sp.pi * t)\n",
    "\n",
    "# 计算 (u·∇)u 的 x 分量\n",
    "term1_ux = ux * sp.diff(ux, x)\n",
    "term2_ux = uy * sp.diff(ux, y)\n",
    "x_component = term1_ux + term2_ux\n",
    "\n",
    "# 计算 (u·∇)u 的 y 分量\n",
    "term1_uy = ux * sp.diff(uy, x)\n",
    "term2_uy = uy * sp.diff(uy, y)\n",
    "y_component = term1_uy + term2_uy\n",
    "\n",
    "# 输出结果\n",
    "# print(\"(u·∇)u 的 x 分量：\")\n",
    "y_component\n",
    "# sp.pprint(x_component)\n",
    "# print(\"\\n(u·∇)u 的 y 分量：\")\n",
    "# sp.pprint(y_component)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Collecting sympy\n",
      "  Downloading sympy-1.14.0-py3-none-any.whl (6.3 MB)\n",
      "\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m6.3/6.3 MB\u001b[0m \u001b[31m4.8 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m00:01\u001b[0m00:01\u001b[0m0m\n",
      "\u001b[?25hCollecting mpmath<1.4,>=1.1.0 (from sympy)\n",
      "  Downloading mpmath-1.3.0-py3-none-any.whl (536 kB)\n",
      "\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m536.2/536.2 kB\u001b[0m \u001b[31m29.2 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
      "\u001b[?25hInstalling collected packages: mpmath, sympy\n",
      "Successfully installed mpmath-1.3.0 sympy-1.14.0\n",
      "Note: you may need to restart the kernel to use updated packages.\n"
     ]
    }
   ],
   "source": [
    "pip install sympy"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
